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Summary. The evaluation of the formative process in the University system has been assum- 
ing an ever increasing importance in the European countries. Within this context the analysis 
of student performance and capabilities plays a fundamental role. In this work we propose a 
multivariate latent growth model for studying the performances of a cohort of students of the 
University of Bologna. The model proposed is innovative since it is composed by: (1) mul- 
tivariate growth models that allow to capture the different dynamics of student performance 
indicators over time and (2) a factor model that allows to measure the general latent student 
capability. The flexibility of the model proposed allows its applications in several fields such as 
socio-economic settings in which personal behaviours are studied by using panel data. 

Keywords: University evaluation, student capability, Data warehouse, longitudinal and mixed 
data, generalized linear latent variable models. 

1. Introduction 

The Bologna Process started in 1999 with the aim of creating a European Higher Education 
Area, in which students could choose from a wide and transparent range of high quality 
courses and benefit from smooth recognition procedures. The Bologna Declaration has ini- 
tiated a series of reforms needed to make European Higher Education more compatible and 
comparable, more competitive and more attractive than before, both for Europeans and for 
students from other continents. Hence, the evaluation of formative processes has received 
a growing attention by policy makers and public agents in order to identify critical factors 
for achievements that can improve curricula, instructional strategies, and conditions for 
learning. 

An important emerging problem is the comparison between students' performances i) when 
different supporting and tutoring actions are adopted during the course of studies, ii) in 
presence of different personal situations. With this purpose, several Universities have cre- 
ated Data WareHouse (DWH) systems to collect detailed multivariate individual responses 
over time, which consist of mixtures of count, categorical, and continuous observations. 
These longitudinal data allow to answer questions about student progress, evaluate how 
each individual performs over time (within-individual change), predict the main differences 
among individuals in their change (interindividual differences in change). However, in pres- 
ence of multidimensional observations, a challenging problem is the characterization of both 
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temporal and cross-sectional dependencies among response variables having different mea- 
surement scales. In such cases, it is natural to consider models in which the dependency 
among the responses is due to the presence of both one or more latent variables and random 
effects, as shown in several approaches developed in the literature. 

Roy and Lin (2000) proposed a 2-step linear mixed model applied to multiple continuous 



outcomes. These authors use time-dependent factors to account for correlations of items 
within time. On the other hand, a random intercept and random effects are introduced for 
explaining the correlations across time of both items and time-dependen t latent variables, 
respectively. An extension to such models is provided bv lDunson (20031 ). who introduced 
a dynamic latent trait model for multidimensional longitudinal data in the context of the 
Generalized Linear Latent Variable Model (GLLVM) so that different kinds of observed 
variables can be considered. An autoregressive structure that allows for covariates is used 
to model the structural part. The model is estimated by using the MCMC procedure. 
Within the same framew ork, a full information likelihood estimation method via the EM 
algorithm is developed bv lCagnone et al~ (2009) with parti cular attention to ordinal data. 
A very general approach is represented by multilevel models (pkrondal and Rabe-Hesketh. 20041 ) 
that allow to deal with longitudinal and/or multidimensional mixed data. In presence of 
repeated measures, occasions are viewed as first level units whereas respondents are second 
level units. With multidimensional data, first level units are represented by items nested 
within individuals. When both the dimensions are considered more complex hierarchical 
structures have to be taken into account. 

Multidimensional and longitudinal data are also treated within the traditional structural 
equation approach (SEM). They are modeled in two diffe rent ways. 

According to the first one (|J5reskog and S5rbom. 20011 ). a standard confirmatory factor 
model is considered and its the main feature is that the corresponding error terms are cor- 
related over time. Moreover, the latent variables are identified by setting to 1 the same 
loading over time. 

The second ap proach is represented by l atent growth models, widely applied in the anal- 
ysis of change (jSinger and Willett. 20031 ). Random effects are included into the model to 
account for the individual differences both in the initial status and in the rate of growth. 
The peculiar feature of such models is that th e random coefficients are treated as latent 
variables within the traditional SEM approach ( Muthen and Khoo. 1998h . In this context, 
univariate analysis are usually performed by studying the temporal dynamics of a single 
indicator, considered as a proxy of the individual performance. Multivariate extensions 
essentially consist of modeling the trajectori es of several items separately, and then t o allow 
for correlations among random coefficients (Boll en and Curran. 20061: iRavkov. 20071) . 
This work is motivated by the study of the data coming from the DWH of the University 
of Bologna. We focus on the achievements of a cohort of students enrolled in 2001 at the 
Faculty of Economics. Multiple items are present in the data set. Their behaviour over time 
can be classically analyzed by means of multivariate latent curves, where all the variability 
between items is captured by the correlation of the random coefficients. However, some of 
these items can be seen as indicators of student latent capabilities. Hence, part of their 
variability can be due to latent constructs. In order to take into account these two impor- 
tant aspects simultaneously we propose a new, general, class of models, that consists of two 
parts: i) multivariate latent curves that describe the behaviour of each item over time, ii) a 
factor model that specifies the relationship between manifest and latent variables. Although 
these two components have been widely developed in the literature separately, the novelty 
of our proposal lies in integrating them into a unique framework. 
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The model is developed within the GLLVM framework. According to this approach, the 
response variables are assumed to follow different distributions of the exponential family, 
with item-specific linear predictors depending on both time-specific covariates, latent vari- 
ables, and measurement errors. Moreover, we extend the GLLVM by including item-specific 
random coefficients so that each item has its own trajectory over time. Our approach has 
clear advantages with respect to both multivariate latent curves since the unobservable ca- 
pabilities of the single individual can be taken into account, and the dynamic factor models 
since we incorporate a more flexible treatment of the temporal dynamics of the items. The 
latter are also assumed to be heteroscedastic over time. 

The paper is organized as follows. In Section 2 we present the data source and perform 
an exploratory analysis to demonstrate the potential of our approach in describing student 
performances over time. Section 3 describes the proposed methodology in terms of model 
specification, identification and estimation. In Section 4 we present the results of model es- 
timation for the overall data set and for different temporal patterns observed in the sample. 
We conclude with a discussion in Section 5. 



2. Data 

The data set analyzed was extracted from the Data warehouse of the University of Bologna. 
This latter is a system that collects and constantly updates informations by integrating data 
coming from sources of different nature. The project started in 2002 in order to support 
planning, control and decision processes. 

The DWH contains a great amount of information per each student and allows to build 
the overall university student career. It is also possible to find socio-demographic informa- 
tion (gender, country/region of origin, etc.) and the mark obtained in the final exam of the 
High School. We decided to analyze the cohort of n = 821 students enrolled at the Faculty 
of Economics in the academic year 2001/2002 since this Faculty is one of the biggest of 
the University of Bologna and such year is the first available in the DWH, so that several 
time points can be observed. Therefore, this data set does not contain missing data in the 




Fig. 1. Trajectories of AM in the three academic years for Overall, Grad, Nograd patterns. 

first three years of the study. After the third year, the presence of missing data is due to 
different reasons: either students that got the degree in time (three years), or drop outs, or 
simply missing information. Hence, we analyzed the performance of the selected students 
in the first three academic years: ti = 2001/2002; t 2 = 2002/2003; t 3 = 2003/2004. In the 
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Table 1. Descriptive statistics of AM in the three academic years 







Overall 






Grad 






Nograd 




ti 


ta 




ti 


ta 


ts 


ti 


ta 


t 3 


Mean 


23.88 


24.14 


23.59 


26.15 


26.33 


26.32 


23.17 


23.46 


22.74 


Std Dev 


4.62 


4.10 


5.26 


1.85 


1.99 


2.05 


4.99 


4.34 


5.65 










Correlations 








ta 


0.41 






0.68 






0.34 






t 3 


0.32 


0.37 




0.67 


0.76 




0.25 


0.29 





Table 2. Number of students per NE over time for the Overall sample, 
Grad and Nograd 
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13 
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3 
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1 
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14 
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data set, two variables are available for evaluating the student performance over time: the 
marks obtained in every exam and the number of exams taken per each time point. The 
average of the marks (AM) per each student in each academic year is considered. In the 
left side of Figure [1] the AM trajectories for the full sample of students are reported. They 
vary from 18 (minimum mark to pass the exam) to 30 (maximum mark) even if they are 
mostly concentrated in the range 21-28, and the overall mean (black line) is over 24 for 
all the observed time points. Within the selected sample we can distinguish two different 
temporal behaviours: the first one concerns students who got the degree in the first three 
years (Grad) while the second one concerns students that at t 3 did not manage to get the 
degree yet (Nograd). Indeed, as shown in the center and in the right side of Figure [U the 
former presents a higher overall average mark and a lower variability than the latter. In 
TablefJJthe descriptive statistics of AM for the overall sample, Grad (ni = 195) and Nograd 
(ri2 = 626) are reported. Grad presents the highest correlations over the three time points. 
On the other hand, Nograd is very similar to the Overall sample in terms of both first 
and second order moments. The number of exams (NE) is a count variable whose range 
is different in the observed time points. Table [2] shows the number of students classified 
according to NE both taken in the three time points and the groups defined before. The 
hyphens indicate that in the first academic year, a student can take at most eight exams. 
We can observe that, in general, Grad students present the same behaviour over the three 
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Table 3. Associations between variables 



Pairs 




Overall 






Grad 






Nograd 






x 2 


df 


p-value 


x 2 


df 


p-value 


x 2 


df 


p-value 


NE tl vs NE t2 


209.10 


30 


0.000 


73.64 


30 


0.000 


80.73 


30 


0.001 


NE tl vs NEt 3 


219.178 


30 


0.000 


85.26 


20 


0.000 


89.51 


30 


0.000 


NE t2 vs NEt 3 


199.02 


36 


0.000 


88.18 


24 


0.000 


110.74 


36 


0.000 


AM tl vs NE tl 


147.45 


15 


0.000 


32.25 


15 


0.006 


90.94 


15 


0.000 


AM; 2 vs NE t2 


92.53 


18 


0.000 


21.39 


18 


0.260 


36.50 


18 


0.006 


AMt 3 vs NEt 3 


137.81 


18 


0.000 


91.52 


9 


0.000 


84.46 


18 


0.001 



years, that is, they take a number of exams greater than three, and concentrate it between 
six and seven. On the contrary, Nograd students take a number of exams equal or greater 
than zero, mostly concentrated between four to six. Moreover, few students of the overall 
sample take more than eight exams at t 2 (only 4,1%) and at t% (only 8,1%). It can be 
useful to evaluate whether i) the variable NE shows a dependence over time and ii) there 
is an association between the variable AM and NE within the same time. To this aim, 
the variable AM has been recoded for all the time points into four classes according to 
the quartiles of the distribution. As for the variable NE, the categories from 1 to 3 and 
categories from 9 to 14 have been collapsed to avoid the problem of sparseness that affects 
these data in the extreme categories. In Table |3] the values of the Chi square tests (with 
associated p-values) are reported for all the pairs of NE over time and for all the pairs of 
AM and NE within time. The association between AM and NE for Grad is not significant 
at time t%, whereas all others are significant, indicating that both the variables can be good 
indicators of the student capability. 

3. Modeling and estimation 

For a given student a record consists of number of exams (NE) and the average marks (AM) 
achieved in each academic year. Hence, in this section we present parametric growth models 
for mixed observations, and provide special treatment to count and continuous responses. 

3. 1. Multivariate latent growth curves 

Suppose that J items are observed for n individuals at Tj,j = 1,...,J, different time 
points. The measured outcomes for a randomly selected individual are denoted by y = 
(yi,:.,yj,—,yj)', where the elements y 5 = (j/y, yrj)', 3 = 1,2,..., J, consist of mix- 
tures of count (j = 1, J\) and continuous (j = J\ + 1,..., J) responses. In analyzing 
data of this type, a challenging problem is the characterization of both the temporal and 
the cross-sectional dependency among variables that have different measurement scales. 
In such cases, it is natural to consider models in which the dependencies are due to the 
presence of both seve ral latent variables and random effects, stacked into the vector rj 
(|Cagnone et al.. 20091) . The marginal distribution of the overall responses is given by 

f(y\x) = J g(y\rj,x)h(ri)dT) (1) 

where g(y\rj, x) is the conditional distribution of the responses y given a set of covariates x 
and the latent variables rj, and h(r)) is their prior density function. We refer to the GLLVM 
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frame work developed in Bartholomew and Knott f!999l ) for multivariate mixed responses 



and in Moustaki (2003) when covariate effects are included. We extend that framework 



to allow for multivariate longitudinal data. One of the main assumption of the GLLVM 
approach is the conditional independence of the responses (within and over time) given the 
latent variables, that is 

g(y\v,x) = YlY[g(y t j\v,x t j) (2) 

3=1 t=l 

where g(ytj\r), x t j) is a distribution of the exponential family. For count data, 

where n t is the number of "trials" (or opportunities for an event). Even if counts are 
generally modeled using a Poisson distribution, when the events being counted for a unit 
occur at a constant rate in continuous time and are mutually independent, a binomial 
distribution is more appropriate. Indeed we will deal with counts corresponding to the 
number of exams that the student takes in each academic year, given the maximum number 
of exams nt observed for that year. The counts have a binomial distribution if the events 
for a unit are independent and equally probable. This assumption is satisfied in our model 
by assuming the conditional independence of the responses given the latent variables rj. 
On the other hand, for continuous data 

g{y tj \ rl ,X tj ) = -^^L\(^i^^\ j = J 1 + l,...,J (4) 

where <jy is the variance of the continuous responses supposed to be heteroscedastic over 
time and between items. 

As in the classical generalized linear model, Vtj is the linear predictor for the jth outcome 
at time t and the link between the linear predictor and the conditional means of the random 
distributions can be any monotonic differentiable function. In this context, the link is the 
logit of the probability associated to each count for the binomial distribution defined in eq. 
([3]) and the identity function in the case of the normal distribution defined in eq. (|4]). 

For both kinds of observed variables the linear predictor is defined as 

p q b 



v t j = ^/3rjA t r + Xkjz k + ^2jtjix t ji t = l,...,Tj, j = 1, J (5) 

k=l 1=1 

V tj = W tj T) + J tj Xtj, j = 1, J (6) 



r=0 k=l 1=1 

and in matrix form 



where 

w ti — [1, At, Af , 0, 0, An, A21, X q i] 



w t j 
7tj 

Xtj 



[0, 0, 1, A t , A t , Ai,/, A 2 j, X q j], 

btji, -,itjb], 

[X t jl, ...,X t jb\. 
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The latent variables r) = (j3oj, /3 p j, Zi, ..., Zq) are random effects and latent traits that 
account for both the temporal and the cross-sectional dependence between items. The 
jtji are fixed regression parameters representing the effects of the covariates Xtji , I — 
1, As it is done with "classical" univariate growth models, the random coefficients 

f3j = (0oj, ....... ,/3 P j), j = 1, J, and the corresponding loadings Xt,t = l,...,Tj, are in- 
troduced in order to describe the temporal behaviour of each item, p being the degree of the 
fitted trajectory. The At's either can be fixed, in the case of linear polynomials, or can be 
parameters to be estimated if a nonlinear function is more appropriate. The model is very 
flexible since it allows to specify different temporal dynamics for each item. The common 
factors z — (zi, . . . , z q ) can represent traits of an individual {e.g. general and specific stu- 
dent abilities, intelligence, etc. ) and determines the correlation between multiple responses 
despite their temporal behaviour. 

By defining the linear predictor in this way, the temporal dependence between items as 
well as the autocorrelation of each item are explained by variance and covariance elements 
related to the random growth parameters f3j,j = 1, J. On the other hand, the cross- 
correlation between items despite their temporal behaviour is caught by the factor model via 
the loadings Xj Z — Xj q ),j = 1, J- These assumptions are contained in the prior 

density function of the latent variables, h(rj), supposed to be a multivariate normal density 
with mean vector /x^ = (^/3 01 , ...,fx 0J)1 , ^/3 0J , H/3 pJ ,0, 0) and covariance matrix = 

*/5 



* 2 



where *b p is the full covariance submatrix related to the random effects 



= 1, •■, J, and is the submatrix of the block matrix related to the latent factors z. 
We assume that the random coefficients (3j and the factors z are independent. Furthermore, 
some constraints must be placed to ensure identifiability of the model based on the observed 
data. In particular, as in classical latent variable models, there is indeterminacy related to 
the scales of the latent factors z ( Joreskog. 1 969). This indeterminacy can be eliminated 



by either setting Xij = 1 or letting the variances of the factors be one, with at least one of 
the loadings constrained to be positive for each factor. 

As a consequence of the structural specification of the model via h(r]), the covariances 
between different responses, seeing in the scale provided by the link function, are given by 

Cov(v t ji,Vt'ji) = W t j0^0w' t , j0 +W jz ^! w' jz 
Cov(v t ji,V t j>i) = Wtj0^0w' tj , +W jz ^ w' jlz 
Cov(vtji,V t >j'i) = w tjf3^pw' t , Jll3 +Wj z ^i 3 w' j , z 



where wtjp and Wj z ,t — 1,2,..., Tj, j — 1, 2, J, indicate the coefficients in Wtj related to 
the random effects — 1, J, and the latent traits z, respectively. 

3.2. Estimation 

The parameters of the model are estimated through the Maximum Likelihood (ML) method 
via the EM algorithm, since the latent variables r) are unobserved. The EM starts with 
initial values of the parameters. The algorithm consists of an expectation and a maximiza- 
tion step. In the expectation step the expected score function from the complete likelihood 
(y,r)\x) given the covariates is computed. In the maximization step updated parameter 
estimates are obtained from the equations derived in the E-step. The whole procedure is 
repeated until convergence. 
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For a random sample of size n, it follows by eq. ([T]) that the complete log-likelihood is 
written as: 

n 

L = ^2^ogf(y i ,r] i \x i ) 

i=l 

n 

= ^2[ l ogg(y i \vi,Xi)+iogh('n i )} (7) 

where g is the likelihood of the data conditional on the covariates, the latent variables and 
the random effects and h is the common distribution function of the latent traits and the 
random effects. From eq. J7]), we see that the first component depends on both the factor 
loadings Xj z ,j = 1, J and the variance parameters er^ = (o\, (Xj,.),j = J\ + 1, J, 
whereas the second component depends on (fx , \&). 



3.2.1. Estimation of y,^ and \f 

From the normality of rj, the second component of the log- likelihood given in eq. |(7J) (up 
to a constant) for an individual i is written as 

log fcfa) = -\ In* - ifa - /i,)* - ^ - /*„)'. (8) 

The expected score function needed for the EM implementation is taken with respect to 
the posterior distribution of the latent variables h(rj i \y i ,Xi). The expected score function 
for the parameter vector \i becomes 

ESi(n v ) = J S i (n v )h(r] l \y l ,x i )dr] i , i=l,...,n (9) 

where 

. dloghOrji) . 
Si(j* n ) = J^ 1 ^ = * (Vi - M„)- 

Similarly, we obtain the score function for that is, 



By solving -^SiOO = and J^iLi ESi(9) = we get explicit solutions for the 

maximum likelihood estimators of fi n and 



3.2.2. Estimation of Xj Z 
The estimation of parameters \j Z ,j 



and <Tj 



1, J and cr^,j = Ji + 1,..., J depends on the 



first component of the log-likelihood given in eq. ([7]). Under the conditional independence 
assumption, the log-likelihood of the count and continuous data can be written as 



logg(yi\vi,Xi 



Ji T 3 

EE 

j=i *=i 



log 


( nt \ 




V Vtji ) 



+ UtjiVtjt - n t log (1 + exp(v t ji)) 



J 



z. E 

jWi + l t=l 



-ilog(27r)-ilog(4)- 



(ytji - vtjif 
24 



(10) 
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The first component refers to count variables and will be used to derive estimates of the 
factor loadings corresponding to such variables, that is, Xj Z ,j — 1,..., Ji- The expected 
score function of the parameter vector Xj z is again taken with respect to the posterior 
HVi\yi,Xi): 



where 



and 



ESi(\j Z ) = J Si(Xj z )h(r] i \y i ,Xi)dr] i , i = l,...,n (11) 
d log g(y t ji\Vi,Xi) 



Si{Xj z ) 



d\ jz 



T 

— -± i * i [ym-rH {1 + eMvtji)) ) , 1 J, (12) 

By replacing eq. (fT2")) into eq. (fTT|) and solving J2™ =1 ESi(\j Z ) — we get non-explicit 
solutions for the parameter vector \j Z . A Newton- Raphson algorithm is used to solve the 
nonlinear maximum likelihood equations. 

From the second component in the likelihood (JTUJ) we estimate factor loadings and vari- 
ance components corresponding to continuous items. The expected score functions for the 
parameters Xj Z ,j = J\ + 1, •••) J, are given by 



ESi(X jz ) = J S l (X JZ )h(q l \y l ,x l )dr] l 



i = 1, n (13) 



where 

Tj t v 

Si(x jz ) = 22 z l^jzi — ^ Zid/tji — 22 PrjXl). 

t = l f=l r=0 

Similarly, we obtain the score function for each element in cr^, j = J\ + 1, J, that is, 



S i(<?tj) = a tj ~ (Vtji ~ V t ji) t = 1, ...,T j 



3 ' 



By solving Yl7=i ESi(Xj Z ) — 0, Y^i=x ESi(cr^) = we get explicit solutions for the maxi- 
mum likelihood estimators of Xj Z and cr'j for j = J\ + 1, J. 

Integrals are approximated by using Gauss-Hermite quadrature points. In order to 
apply the Gauss-Hermite approximation to the integral of equations ©, (JTTJ, and (fT^|) we 
consider t he Cholesky decompo sition of the covariance matrix given by \I> = CC' . As 
shown by ICagnone et al. (2009), this is necessary because the non null submatrices of \&, 



namely and \& z , are not diagonal. 

The steps of the EM algorithm are defined as follows: 

Step 1: choose initial estimates for the model parameters. Starting values for the loadings 
are obtained by fitting separate confirmatory factor analysis models at each time 
points. Initial values for the other parameters are chosen arbitrarily. 

Step 2: compute the expected score functions for all the parameters (E-step). 

Step 3: obtain improved estimates for the parameters by solving the nonlinear maximum 
likelihood equations for the parameters corresponding to the count items and explicit 
solutions for the parameters of the continuous items and the latent distribution (M 
step). 

Step 4: repeat steps 2-3 until convergence is attained. 
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4. Results 

We start the analysis by estimating a model for the overall dataset of students observed 
at the three different time points. As already discussed, the aim of the analysis is twofold: 
(1) analyze the student careers over time with respect to the Number of Exams taken 
in each occasion (NE) and the corresponding Average Marks (AM), and (2) measure the 
general latent capability of the students. Therefore, we analyze how the variables NE 
and AM change over time by means of multivariate latent growth models, and we extend 
these models by including a common factor that can explain the atemporal variability that 
exists between the two items. In particular, since only three different academic years are 
considered, a linear polynomial model (p = 1) could be appropriate to describe the temporal 
pattern of both NE and AM. Measurement invariance over time of the loadings in the one- 
factor model (q = 1) is assumed. Thus, the estimated model (denoted as Model A) is 
characterized by the following linear predictor 

v tj = /3 0j + {t - 1) + XjZ = w tj r], t = 1, 2, 3, j = NE, AM 

where 

w tNE = [1,0,1,2,0,0,0,0^^] 
wtAM = [0,0,0,0,1,0,1,2,A A m], 

V = [Pone, Pine , Pa am , Pi am , A ■ 

Binomial-logistic and Normal heteroscedastic linear regression models are estimated for 
the NE and AM, respectively. The multivariate normal density of the latent variables 
h(rj) has mean vector n v = (m^ojvbjM^ijvb)M/3damjM/3iam>°) and covariance matrix * = 
/ \ff g o \ 

I p l ) ' ^ or identification reasons, the variance of the common latent factor z is set 
equal to 1. 

FORTRAN and R codes have been implemented to estimate the model. (They are available 
from the authors upon request). Parameter estimates of Model A for the overall dataset 
are reported in Table [U 

It can be noticed that both NE and AM present similar temporal dynamics, even if ex- 
pressed on different scales. In terms of the population mean trajectory, NE presents a mean 
initial status equal to 0.249, indicating that in the first year the students take, in mean, 
around 4.5 exams, as expressed in the original scale. The students' progress is described by 
the slope mean parameter (Jtpi, equal to -0.443, reflecting the linear, term-by-term, worsen- 
ing in mean achievement during the second and third years. Both the mean initial status 
and rate of growth are coherent with the descriptive analyses we performed in Section 2, 
since on average the number of exams, as derived by Table [H were 4.76, 4.43, and 4.47 
in ti,t2,t3, respectively. Similarly, at the initial status, students obtain an average mark 
(in mean) around 23.98, but this mean worsens over time as indicated by the mean slope 
parameter jlfj l AM equal to —0.133. 

By looking at the covariances specific of each item in ^ p, students present a higher variabil- 
ity in the initial status than in the rate of growth, with a negative correlation between initial 
status and slope, for both NE and AM. Multivariate latent curves also allow to analyze the 
covariation between the temporal dynamics of NE and AM by estimating cross-covariances 
between random intercepts and slopes of the two curves. There are positive and significant 
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Table 4. Estimates for the overall data set (standard errors in brackets) 



Model A 

Coefficients 



Estimates 



Multivariate growth model 

fc/3 AM 
*3 



Factor model 

Xne 

Xam 

2.2 



a lAM 
a 2AM 
"SAM 



0.249 (0.047) 
0.443 (0.026) 
23.98 (0.326) 
0.113 (0.244) 
/ 0.231* 

-0.107* 0.093* 
1.004* -0.531* 
-0.363* 0.413* 



4.980* 
-2.183* 



2.164* / 



nificant at 5% lc 



0.524 (0.073) 
2.581 (0.449) 

10.080 (0.340) 
9.953 (0.176) 

14.568 (0.395) 



covariances between Pone and /3oam as well as between the random slopes, (3ine and (3\am, 
indicating that students with higher (smaller) average marks in the first year tend to take 
a higher (smaller) number of exams at t\ , and that students with positive (negative) slopes 
for AM generally present a similar pattern for NE. On the other hand, negative covariances 
are estimated between Pone and Piam as well as between Poam and Pine- Differently 
from classical multivariate latent growth modeling, these cross-covariances between NE and 
AM curves are free from the effect of a common latent factor z we estimated via integrating 
the growth curves with a one factor model. 

When manifest variables are of different types, care is needed in the interpretation of 
the factor loadings, depending on the scale of the ytj- In order to interpret the latent 
factor z we shall therefore have to ensure that the As are calibrated so that they may be 
meaningfully compared across variable types. This m ay b e done in a variety of ways but w e 
follow the approach of Takane and De Leeuw (1987 ) and Bartholomew and Knott (19991 ) . 



which provides a parametrization that keeps the interpretation as close as possible to the 
familiar methods of traditional factor analysis. This approach is based on a standardization 
of the coefficients of the latent variable z in order to express correlation coefficients between 
the manifest variable y t j and the factor z. 

For the normal item, Xam denotes the covariance between the manifest variables ytAM and 
the factor z. By dividing Xam by the square root of the variance of the continuous variable 
VtAM, t = 1, 2, 3, we obtain the correlation between ytAM and z, that is 



A tAM — I = = 1—1,4,6. 

+ ^ + (* " 1) 2 ^ 1AM + Ai 



'AM ~ u tAM 



(14) 
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Notice that the correlation varies over time, hence 





= 0.554 


*2AM 


= 0.586 


^3AM 


= 0.505 



The amplitude of the factor loadings is quite similar in all the three occasions and on average 
is equal to 0.548, indicating that the measurement invariance assumption is appropriate. 
On the other hand, for the binomial item , the standardization follows that proposed for 



binary items (jMoustaki and Knott. 20011 ) a nd based on the equivalence of the response 



function and underlying variable approaches (Ta kane and De Leeuw. 19871) . In this context, 



the correlation between a normal variable supposed to be underlying the binomial discrete 
observations and the latent variable z is given by 

Kne= i t= 1,2,3. (15) 

V^/w + 2 ^ - Wfh**,*** + (* " 1 ) 2 ^} 1NE + x %e + 1 

The estimated standardized binomial loadings are 



A 1NE 


= 0.427 


^2NE 


= 0.445 


^3NE 


= 0.435 



which are really close to each other; the amplitude of the standardized factor loadings is on 
average 0.436. The standardized coefficients given for normal and binomial variables can be 
used for a unified interpretation of the loadings, bringing the interpretation close to factor 
analysis. The common factor z explains the interrelationships between the two observed 
items net from their temporal dependence. Both variables are significant indicators of 
this latent capability. Moreover, they both influence positively this unobserved construct. 
Ignoring the presence of a common factor in multivariate latent growth models can lead to 
an overestimation of the cross-variation among multiple curves. 

T he goodness of fit of Model A has been checked separately for the count and continuous 



part (jMoustaki and Knott. 20011 ). As for the count part of the model, significant informa- 
tion concerning the goodness of fit can be found in the margins. In particular, the one-way 
margins of the differences between the observed (O) and expected (E) frequencies under 
the model are investigated; any large discrepancies will suggest that the model does not 
fit well these counts. The Chi square t est as well as high- way margins are not appropriate 



because of the sparseness of the data (jReiser. 19961 ) . Table [5l gives the GF-fit measures, 
calculated as ((O — E) 2 /E), for each Binomial variable ([Bartholomew et al.. 20021 ). We can 



observe that the GF-fits are not good, especially those on count 7 for Uine, on 8 and 9 for 
V2NE, and on counts 6 and 7 for y^NE- Reasons of this misfitting of Model A on the overall 
sample will be next investigated. 

For the normal part of the model we check the discrepancies between the sample correla- 
tion matrix and the one estimated from the model, as illustrated in Table[6]for the variables 
DiAMi DiAM) an d U3AM- The discrepancies between observed correlations and those esti- 
mated are particularly small, indicating that the fit of the model for the normal variables 
is good. 
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Table 5. Count items: GF-fit values 
for the one-way margins. 



Counts 


yiEN 




7/Q R 1 AT 


o 


11.58 


11.65 


11.23 


1 


0.10 


1.59 


3.05 


2 


5.78 


0.22 


10.39 


3 


5.45 


0.39 


15.20 


4 


8.73 


4.97 


5.91 


5 


1.10 


16.18 


0.95 


6 


7.58 


5.12 


34.35 


7 


24.47 


0.33 


37.70 


8 


3.65 


23.75 


2.33 


9 




31.37 


0.04 


10 




13.06 


0.66 


11 




0.31 


7.94 


12 




0.02 





13 






0.41 


14 






1.23 



Table 6. Normal items: discrep- 
ancies between sample and esti- 
mated correlation matrices. 





VlEN 


J/2BJV 


VZEN 


yiEN 


0.00 


-0.05 


0.02 


V2EN 


-0.05 


0.00 


-0.05 


V3EN 


0.02 


-0.05 


0.00 



In Section 2, we showed that among the 821 students, two different temporal patterns 
were evident, one related to those students who graduated at £3 (Grad) and the other to 
students who did not get the degree regularly at the third year (Nograd). Hence, we shall 
analyze these two different groups of students in order to investigate the reasons of the 
poorness of fit for the count part of the Model A in the overall dataset. Therefore, in the 
following, we fit Model A to the Grad and Nograd students, separately. 

4. 1. The Graduate students 

We first consider the 195 students who graduated at t 3 , and we start by fitting the Model 
A described above. The results of the estimation are reported in Table [7) It can be noticed 
that parameter estimates corresponding to both the multivariate latent growth and the 
factor parts of Model A differ substantially from what we obtained for the overall sample. 
As for the former, both NE and AM present higher mean initial status than the overall 
sample. They are equal to 6.05 for NE and to 26.17 for AM, as expressed in the original 
scales. Furthermore, the mean trajectory for NE has a worsening pattern over time, whereas 
the one corresponding to AM presents an increasing but not significant temporal behaviour. 
By looking at the variability around the mean trajectories, this is significant in the initial 
status (0.133) and in the rate of growth (0.072) corresponding to NE which also shows a 
negative correlation between the random intercept and slope. On the other hand, there is 
not a significant variability for AM with respect to its mean trajectory, indicating that the 
Grad students show a similar pattern over time. This finding is also evident in all the (not 
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significant) covariances between random coefficients of NE and AM. 

As for the factor part, the loading associated to the binomial variable is very close to 
and not significant, suggesting that the number of exams for these students is not a 
measure of their capability. However, also the loading associated to AM is not significant. 
This means that for Grad students it makes no sense to specify a common factor related to 
the variables AM and NE. A justification could be found from the not significant association 
between AM and NE for Grad in time t2, as shown in Section 2. This evidence can be a 
hint to test the assumption of measurement invariance over time. If we estimate a model 
where this assumption is relaxed, denote it as Model B (Table [7]), we can notice that the 
time dependent loadings are very different. This is particularly true for the variable NE, for 
which the loading does not change greatly in the first two time points but becomes negative 
in £3. For the variable AM the loading increases over time. Clearly, the measurement 
invariance cannot be assumed. This result is also confirmed by the AIC and BIC criteria 
that show that Model B is better than Model A. However, from our viewpoint, Model B is 
meaningless in the factor part. Moreover, if we look at the one-way margins associated to 
Model B (Table E|) we can see that again there are goodness of fit problems at time points 
and £3. 

These results for Grad highlight two different aspects of the analysis. First of all, we have 
a slight individual variability around NE and AM mean trajectories, indicating a similar 
temporal behaviour of these students. Moreover, the higher values of these means compared 
to those obtained for the overall sample show a good performance of this group. Secondly, 
in this case the two variables are not measures of a latent construct that in the overall 
sample we identified as capability. 



4.2. The Nograduate students 

Coherently with the previous analysis, we first estimated Model A for Nograd students. 
The results are reported in Table [HI The growth model shows results similar to the overall 
sample with a higher variability for AM with respect to its mean trajectory. In this case the 
results related to the factor model are very interesting. Differently from what we found for 
the Grad group, the loadings are both significant and positively related to the latent variable 
and indicate that, as in the overall data set, the factor model is appropriate. Differently 
from what we found in all previous analysis, the GF-fits of this model are satisfactory for 
all the observed time points, as reported in Table [TUl Thus the count part of the model 
is well fitted by the binomial distribution. Also for the normal part the fit is very good, 
the discrepancies between observed and estimated correlations being very low (Table fTTj) . 
Therefore, Model A fits well the Nograd students data. 

If we look again at Table [9] Model A, it can be noticed that the values of <j\ m are quite 
similar over time; thus, it can be interesting to evaluate if AM is homoscedastic over time. 
The results of the estimation of the model with homoscedastic errors (Model C) are reported 
in the bottom part of Table[9] Although all the parameter estimates do not change abruptly, 
the AIC and BIC are slightly better for Model A than for Model C, suggesting that such 
assumption does not hold. Thus, the comparison between the loadings of AM and NE for 
evaluating the influence of each item on the latent variable requires their standardization 
according to eq. (|14[) for AM and to eq. (fT5)) for NE. We get the following standardized 
loadings: \ AM = (0.525,0.560,0.534), and \ NE = (0.401,0.414,0.388). As in the overall 
data set, the correlation between z and AM is slightly higher than that between z and NE. 
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Table 7. Estimates for graduate (Grad) students 



Model A 

Coefficients 



Estimates 



Multivariate growth model 

~ 2 

^PqNE 

~2 

l^pxNE 
'2 

~ 2 

M/3i AM 



Factor model 
A 



NE 

Aam 

-2 



AIC= 4304.090 
BIC= 4309.602 



1.133 (0.086) 
-0.577 (0.050) 
26.172 (0.149) 
0.080 (0.061) 



( 0.133* 

-0.097* 0.072* 

0.441 -0.327 

\ 0.044 -0.026 

*: significant at 5% level. 



0.005 (0.262) 
0.670 (3.360) 
0.702 (0.147) 
1.292 (0.113) 
0.328 (0.185) 



2.197 
-0.086 



0.376 / 



Model B 

Coefficients 



Estimates 



Multivariate growth model 

~ 2 

^BnNE 
-2 

H-fixNE 
-2 

H-/3 AM 
~2 

f'pxAM 



Factor model 
\ine 

^2NE 

A3JVB 

AlAM 

\2AM 

\ZAM 
^2 

"1AM 
~ 2 

G 2AM 
°~3AM 



1.193 (0.093) 
-0.594 (0.066) 
26.225 (0.192) 
0.188 (0.079) 



/ 0.136* 

-0.103* 0.078* 

0.445 -0.340 

0.025 -0.021 

significant at 5% level. 



2.727 
-0.167 



0.317 / 



0.300 (0.413) 
0.465 (0.594) 
-0.103 (0.185) 
0.156 (0.339) 
0.440 (0.571) 
0.734 (0.886) 
0.616 (0.129) 
1.300 (0.114) 
0.351 (0.177) 



AIC= 4226.847 
BIC= 4273.520 
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Table 8. Count items: GF-fit values 
for the one-way margins, Model B. 



Counts 


yiEN 


J/2BJV 


VSEN 


2 


0.34 


0.05 




3 


4.65 






4 


7.46 


0.03 




5 


4.20 


16.28 


18.16 


6 


1.37 


58.32 


11.10 


7 


5.29 


13.48 


19.52 


8 


0.86 


15.12 


0.00 


9 




23.81 


0.40 


10 




14.41 


0.50 


11 




0.02 


2.35 


12 




1.38- 




13 








14 






30.51 



5. Discussion 

In this paper we extended and applied multivariate latent growth models to the analysis of 
student record data collected repeatedly in the Data warehouse system of the University of 
Bologna. The proposed approach is innovative since it allows to evaluate both the student 
performance over time and individual capabilities simultaneously. Key features include i) 
a flexible modeling of the temporal dynamics of the observed variables via specific latent 
curves, and ii) an extension of the multivariate growth model that incorporates a factor 
part. Such component explains the association between the observed items by means of 
latent variables, interpreted as different traits or capabilities. 

The complexity of the model proposed lies in different aspects, such as the presence of mixed 
data, the possibility of both including several latent variables/random effects and estimating 
specific temporal patterns for the observed variables. Hence, computational problems occur 
in the parameter estimation. We successfully solved them by implementing an ad hoc EM 
algorithm (Fortran and R code are available upon request by the authors). As far as we 
know, commercial software does not allow to treat all these aspects simultaneously. 

We demonstrated, via different specifications of the model, how our general approach 
can provide insights into the data structure. In particular, the analysis carried out on a 
cohort of students enrolled at the Faculty of Economics observed at three different time 
points highlighted an heterogeneity in the overall data set in terms of both average marks 
and number of exams. This is due to the presence of different temporal patterns within the 
cohort, since we have students who regularly graduate at £3 (Grad) and students who did 
not manage to get the degree within the third year (Nograd). Grad students perform very 
well in terms of both number of exams and average marks with similar temporal pattern. 
Nograd take a lower number of exams with lower average marks, but within this group we 
have a significant variability both in the initial status and in the rate of growth. The factor 
part of the model for the Grad student fails in measuring a general capability by means of 
the observed indicators considered. We found that this fact depends on the fundamental 
assumption of measurement invariance of items over time. Such assumption does not hold 
in this case. On the contrary, the model fits well the data of Nograd students. What we 
called atemporal latent capability is well measured by the average mark and the number of 
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Table 9. Estimates for Undergraduate (Nograd) students 



Model A 

Coefficients 



Estimates 



Multivariate growth model 

~2 

t L Q NE 

~2 

L l 1 NE 

~2 

^0 o AM 

~2 

A* 0i AM 



Factor model 

~2 

"1AM 
a 2AM 
°~3AM 

AIC= 18773.348 
BIC= 18788.483 



-0.030 (0.055) 
-0.353 (0.035) 
23.337 (0.431) 
-0.183 (0.326) 

/ 0.171* 

-0.109* 0.131* 

0.930* -0.610* 

\ -0.552* 0.658* 

*: significant at 5% level. 



0.473 (0.119) 
2.609 (0.706) 
12.356 (0.567) 
12.484 (0.269) 
17.307 (0.612) 



5.514* 
-3.304* 



3.492* / 



Model C 

Coefficients 



Estimates 



Multivariate growth model 

'2 

t L o NE 

"2 

L l 1 NE 

"2 

1^00 AM 

"2 

A* 0! AM 



-0.031 (0.051) 
-0.353 (0.032) 
23.346 (0.406) 
-0.215 (0.294) 



/ 0.184* 
-0.117* 
0.935* 
\ -0.570* 



0.133* 

-0.619* 5.180* 
0.709* -3.241* 4.052* 



*: significant at 5% level. 



Factor model 
Aam 

"2 

° AM 



0.464 (0.079) 
2.502 (0.620) 
13.642 (0.234) 



AIC= 18777.719 
BIC= 18791.260 
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Table 10. Count items: GF-fit val- 
ues for the one-way margins, Model A, 
Nograd. 



Counts 


yiEN 


1J2EN 


y3EN 





7.69 


11.56 


5.01 


1 


0.86 


1.67 


4.91 


2 


7.91 


0.18 


9.11 


3 


3.48 


1.33 


8.28 


4 


2.16 


7.46 


0.17 


5 


1.01 


12.83 


2.96 


6 


9.86 


0.33 


20.67 


7 


12.61 


7.89 


8.47 


8 


0.27 


10.92 


0.41 


9 




14.13 


0.27 


10 




5.54 


3.57 


11 






3.21 


12 




0.35 




13 






1.53 


Table 11. 


Normal items: 


discrep- 


ancies between sample and esti- 


mated correlation matrices, Model 


A, Nograd 












yiEN 


y'2EN 


yZEN 




Vien 


0.00 


-0.07 


0.02 




U2EN 


-0.07 


0.00 


-0.08 




l/SEN 


0.02 


-0.08 


0.00 



exams taken, both being significantly related to the latent variable. The good performance 
of the model is confirmed by the analysis of some goodness of fit statistics. 
The heterogeneity observed in the patterns of graduate and nograduate students has im- 
plications on the results of the overall data set. On the one hand parameter estimates for 
the overall sample are quite similar to those of Nograd students, the factor loadings of the 
two variables being both significant. This is in part due to the larger sample size of Nograd 
students. However, the presence of different performances of the Grad students reflects on 
the poorness of fit of the count variable in the overall data set. 

The model proposed here was motivated by the study of the students' achievements and 
the good results obtained clearly show its appropriateness. However, such methodology 
can be applied successfully in many other fields, such as socio-economic settings in which 
personal behaviours are studied by using panel data collected through the administration 
of questionnaires. 

In the present example no covariates have been considered. In practice, we may have 
useful time dependent and time independent covariates such as gender, region of origin, age, 
etc. that can be incorporated into the model. In particular, an emerging field of investiga- 
tion is based on the comparison of th e performances of students who completed the degre e 
compared with those who abandoned (jSmith and Navlor. 2001tlDraper and Gittoes. 2 004). 
Furthermore, the treatment of missing data would allow to extend our analysis to more time 
points and evaluate if non linear or higher degree polynomial trajectories can describe the 
temporal behaviour of the items s tudied. Preliminary studies performed with the software 
LISREL (jBianconcini et. al. 20071 ) showed how different latent curves can fit the weighted 
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average marks for different groups of students within the cohort analyzed. Such problems 
will motivate our future investigations along these lines of research. 
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